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ABSTRACT 

The excursion set theory based on spherical or ehipsoidal gravitational collapse pro- 
vides an elegant analytic framework for calculating the mass function and the large- 
scale bias of dark matter haloes. This theory assumes that the perturbed density field 
evolves stochastically with the smoothing scale and exhibits Markovian random walks 
in the presence of a density barrier. Here we derive an analytic expression for the 
halo bias in a new theoretical model that incorporates non-Markovian extension of 
the excursion set theory with a stochastic barrier. This model allows us to handle 
non-Markovian random walks and to calculate perturbativly these corrections to the 
standard Markovian predictions for the halo mass function and halo bias. Our model 
contains only two parameters: k, which parameterizes the degree of non-Markovianity 
and whose exact value depends on the shape of the filter function used to smooth 
the density field, and a, which parameterizes the degree of stochasticity of the barrier. 
Appropriate choices of k and a in our new model can lead to a closer match to both the 
halo mass function and halo bias in the latest A^-body simulations than the standard 
excursion set theory. 
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Dark matter haloes typically form at sites of high density 
peaks. The spatial distribution of dark matter haloes is 
therefore a biased tracer of the underlying mass distribution. 
A standard way to quantify this difference between haloes 
and mass is to use a halo bias parameter bh, which can be 
defined as the ratio of the overdensity of haloes to mass, or 
as the square root of the ratio of the two-point correlation 
function (or power spectrum) of haloes to mass. 

Like the halo mass function, analytic expressions for 
the halo bias can be obtained from the excursion set the- 
ory (Bond et al. 1991) based on the spherical gravitational 
collapse model (Cole & Kaiser 1989; Mo & White 1996). In 
the excursion set theory, the density perturbation evolves 
stochastically with the smoothing scale, and the problem of 
computing the probability of halo formation is mapped into 
the first-passage time problem in the presence of a (constant) 
barrier. The approach to the clustering evolution is based 
on a generalization of the peak-background split scheme 
(Bardeen et al. 1986), which basically consists in splitting 



the mass perturbations into a fine-grained (peak) compo- 
nent filtered on a scale R and a coarse-grained (background) 
component filtered on a scale _Ro 3> R- The underlying idea 
is to ascribe the collapse of objects on small scales to the 
high frequency modes of the density fields, while the action 
of large-scale structures of these non-linear condensations is 
due to a shift of the local background density. 

Comparison with A'^-body simulations finds that the 
spherical collapse model underpredicts the halo bias for low 
mass halos (Jing 1998; Sheth & Tormen 1999). The discrep- 
ancy reaches a factor of ~ 2 at Af ~ O.OIA/* , where Af» is the 
characteristic nonlinear mass scale (defined by cr(M«) — 1 
where a'^{M) is the variance of the density field in a vol- 
ume of radius R containing the mass M). Sheth et al (2001) 
obtained an improved formula for the halo bias by using 
a moving barrier whose scale-dependent shape is motivated 
by the ellipsoidal gravitational collapse model. Compared to 
the spherical collapse model, this formula predicts a lower 
bias at the high mass end and a higher bias at the low mass 
end (see Fig. 1 below). The resulting bias is shown to be too 
high at the low mass end by ~ 20% compared with simu- 
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lation results. Further modifications have been introduced 
that either used the functional form of Sheth & Tormen 
(1999) or Sheth et al (2001) with new fitting parameters 
(e.g., Tinker et al 2005), or proposed new fitting forms al- 
together (e.g., Seljak & Warren 2004; Pillepich et al 2010; 
Tinker et al 2010). 

Our goal in this paper is not to improve on the accu- 
racy of the fits to the halo bias, but rather to gain deeper 
theoretical insight by deriving an analytical expression for 
the halo bias using a new model. This model modifies the 
excursion set theory by incorporating non-Markovian ran- 
dom walks in the presence of a stochastic barrier. It is based 
on a path integral formulation introduced in Maggiore & 
Riotto (2010a, b), which provides an analytic framework for 
calculating perturbatively the non-Markovian corrections to 
the standard version of the excursion set theory. Mathemat- 
ically, the non-Markovianity in the theory is related to the 
choice of the filter function necessary to smooth out the 
density contrast. As soon as the filter function is different 
from a step (tophat) function in momentum space, the ex- 
cursion of the smoothed density contrast is non-Markovian, 
namely, every step depends on the previous ones and the 
random walk acquires memory. As the computation of the 
bias parameter hh amounts to computing the first crossing 
rate with a non-trivial initial condition at a large, but not 
infinite, radius, the non-Markovianity makes the calculation 
much harder than the Markovian case. 

Furthermore, the critical value for collapse in our model 
is itself assumed to be a stochastic variable, whose scatter 
reflects a number of complicated aspects of the underlying 
dynamics. The gravitational collapse of haloes is a complex 
dynamical phenomenon, and modeling it as spherical, or 
even as ellipsoidal, is a significant oversimplification. In ad- 
dition, the very definition of what is a dark matter halo, both 
in simulations and observationally, is a non-trivial problem. 
Maggiore & Riotto (2010b) proposed that some of the phys- 
ical complications inherent to a realistic description of halo 
formation can be included in the excursion set theory frame- 
work, at least at an effective level, by taking into account 
that the critical value for collapse is itself a stochastic vari- 
able. 

In Section 2 we review briefiy the derivation for the 
halo mass function in the Markovian excursion set theory 
(Sec 2.1) and the path integral approach used to introduce 
non-Markovian terms (Sec 2.2) and stochastic barriers (Sec 
2.3) into the theory. Section 3 is devoted to the discussion 
of halo bias, including a review of the standard derivation 
in the Markovian case (Sec 3.1), and a summary of our new 
derivation in the non-Markovian model (Sec 3.2 and 3.3). 
The details of how the halo bias is calculated from the condi- 
tional probability for two barrier crossings in the new model 
is provided in the Appendix. In Section 4 we compare the 
predictions for the halo bias and mass function in our new 
model with those from the Markovian model (both spherical 
and ellipsoidal collapse) and A^'-body simulations. 



troduce non-Markovian terms in the presence of a stochastic 
barrier (Maggiore & Riotto 2010a, b). 



2.1 Brief review of the excursion set theory 

The basic variable is the smoothed density contrast, 

S{yi,R)= / dViy(lx-x'l,i?)5(x'), 



(1) 



where 5(x) = p{'x.)/p — 1 is the density contrast about the 
mean mass density p of the universe, W^(|x — x'|, R) is the 
filter function, and R is the smoothing scale. We are inter- 
ested in the evolution of (5(x, R) with smoothing scale i? at a 
fixed point x in space, so we suppress the argument x from 
this point on. It is also convenient to use, instead of R, the 
variance S of the smoothed density field defined by 



S{R) = a{R) 






P{k)W\k,R), 



(2) 



where P{k) is the power spectrum of the matter density 
fluctuations in the cosmological model under consideration, 
and W is the Fourier transform of the filter function W . 
A smoothing radius R = oo corresponds to S = and, 
in hierarchical models of structure formation such as the 
ACDM model, S is a monotonically decreasing function of 
R. We can therefore use S and R interchangeably and denote 
our basic variable by 5{S). 

If the filter function is taken to be a tophat in momen- 
tum space, S{S) then satisfies a simple Langevin equation 
with S playing the role of a "pseudo-time" (Bond et al. 1991) 



dSjS) 
dS 



^ vis) , 



(3) 



where rj{S) represents a stochastic "pseudo-force" whose 
two-point correlation statistic obeys a Dirac-delta function: 



{viSi)niS2)) = 5D{Si-S2). 



(4) 



It then follows that the function II((5o; 5; S), which gives the 
probability density of reaching a value S at "time" S starting 
at a value 5o at S — 0, satisfles the Fokker-Planck equation 



9n _ 1 d^ 
'dS ^ 2~dP 



(5) 



In Bond et al. (1991) this equation was supplemented by the 
boundary condition 



n(5, 5)1,^,^=0 



(6) 



to eliminate the trajectories that have reached the critical 
value 5c for collapse. 

The corresponding solution of the Fokker-Planck equa- 
tion is 



n(5o;<5;S) 



g-(i-io)V(2S) _ ^-{2S,-So-Sy/(2S) 



•(7) 



The probability J^{S)dS of first crossing the threshold den- 
sity Sc between "time" S and S + dS is then given by 



2 NON-MARKOVIAN EXTENSION AND 
STOCHASTIC BARRIER 

In this section we review the main points of the excursion 
set theory (Bond et al. 1991) and then summarize how to in- 



ns) = 



dS 



m 

ds 



2^S3/2 



-S^/(2S) 



(8) 



where we have set (5o = 0. The number density of virialized 
objects with mass between M and M + dM is related to the 
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first crossing probability between S and S + dS by 



dnjM) 
dM 



dM 



jjHS)dS 



2 5, 



£-^2/(2^2) J^ 



dlntj 



M2 dlnM 



■dM 



/M 



diner 



M2 dlnM 



-dM , 



(9) 



where a — S^' is defined in equation (2) and u = 5c/cr. 
Tills expression reproduces the mass function of Press & 
Schechter (1974), including the correct overall normaliza- 
tion that had to be adjusted by hand in Press & Schechter 
(1974). We will use the dimensionless function fiy) defined 
in equation (9) to denote the halo mass function below. 

2.2 Non-Markovian extension 

As already discussed in Bond et al. (1991), a difficulty of 
the excursion set approach in Sec 2.1 is that an unambigu- 
ous relation between the smoothing radius R and the mass 
M of the corresponding collapsed halo only exists when the 
filter function is a tophat in coordinate space: M{R) = 
(4/3)7r iJ'^p. For all other filter functions (e.g., tophat in mo- 
mentum space, Gaussian), it is impossible to associate a 
well-defined mass M{R) (see also the recent review Zentner 
2007). More importantly, S{S) obeys a Langevin equation 
with a Dirac delta noise as in eqs. (3) and (4) only when the 
filter function is a tophat in momentum space. Otherwise, 
the evolution of 5 with the smoothing scale becomes non- 
Markovian, and the distribution function Il{So; S; S) of the 
trajectories no longer obeys the Fokker-Planck equation, nor 
any local generalization of it. In this case, Il{5o; S; S) obeys 
a complicated equation that is non-local with respect to the 
variable S (Maggiore & Riotto 2010a). 

To deal with this problem, Maggiore & Riotto (2010a) 
proposed a "microscopic" approach, in which one computes 
the probability associated with each trajectory S{S), and 
sums over all relevant trajectories. As with any path integral 
formulation, it is convenient to discretize the time variable 
and to take the continuum limit at the end. Therefore we 
discretize the interval [0, S] in steps AS = e, so Sk ~ ke 
with fc = 1, . . . n, and Sn = S, and a trajectory is defined by 
the collection of values {5i, . . . , Sn}, such that 5{Sk) ~ Sk- 
All trajectories start at a value So at time S — 0. 

The basic quantity in this approach is the probability 
density in the space of trajectories, defined as 



W{So;Si 



, Un'i On 



(fc(5(5i)-5i) . . . Sd{S{S„)-S„)){10) 



where 5d denotes the Dirac delta function. In terms of W 
we define 

n4So;S„;Sn)= f dSi .. .dS„-iW (So; Si, ... ,S„;S„) (11) 

where Sn ~ ne, and IIe((5o; S; S) is the probability density of 
arriving at the "position" 5 in a "time" S, starting from Sq 
at time So = 0, through trajectories that never exceeded Sc. 
The problem of computing the distribution function of ex- 
cursion set theory is therefore mapped into the computation 
of a path integral with a boundary at S = Sc. 

The probability density W can be computed in terms 
of the connected correlators of the theory. When the density 



field (5 is a Gaussian random variable, only the two-point 
connected function is non-zero, and one finds 



W{So;Si,...,S„;S„)^ 

dXi d\„ ^y^" >ii'5i-iy^" Xi\i{SiSi)c 



(12) 



27r 



27r 



and Si = S{Si). We will restrict the discussion here to the 
Gaussian case since higher-order connected correlators must 
be included in the non-Gaussian case (Maggiore & Riotto 
2010c,d). 

Consider first the case of a tophat filter in momen- 
tum space, so the evolution of S{S) is Markovian and obeys 
eqs. (3) and (4). Then one can show that the connected 
two-point correlator is given by 



{S{S,)S{Sj))c^ minis,, Sj), 



(13) 



and the integrals over dXi, . . . , d\n in eq. (12) can be per- 
formed explicitly to give 



W'"^{So;Si 



, On ; On) 



=-* e::o'(^-+i-*')' 



(14) 



(27re)"/2 

where the superscript "gm" stands for "Gaussian and 
Markovian." Inserting this expression into eq. (11) it can be 
shown (Maggiore & Riotto 2010a) that, in the continuum 
limit, the corresponding distribution function nf™Q((5o; 5; 5) 
satisfies the Fokker-Planck equation (5) as well as the 
boundary condition (6), and therefore we recover the stan- 
dard result (7) of excursion set theory. 

The interesting case is to generalize the computation 
above to filter functions different from the conventional 
tophat in momentum space. The two-point correlator de- 
pends on the filter function. For a Gaussian filter and a 
tophat filter in coordinate space, for instance, we find 



{S{S,)S{Sj)) = min(S'„ S,) + A{S„ Sj 



(15) 



where A{Si,Sj) — A{Sj,Si) and, for Si ^ Sj, the function 
A{Si, Sj) is well approximated by 



A{Si,Sj) ~ K 






(16) 



with K ~ 0.35 for a Gaussian filter and k ~ 0.44 for a tophat 
filter in coordinate space. The parameter k gives a measure 
of the non-Markovianity of the stochastic process, and the 
computation of the distribution function Il{So; S; S) can be 
performed order by order in n. The technique necessary for 
evaluating the path integral in eq. (11) to first order in k has 
been developed in Maggiore & Riotto (2010a), and will be 
further discussed below. To first order in the non-Markovian 
corrections, the resulting first-crossing rate becomes 



ns) 



1 



2^ 53/2 



^-Si/{2S) _|_ 



2V2^S'^/- 



-r 



25 



(17) 



where r(0, z) is the incomplete Gamma function. For k = 
one recovers the Markovian result in equation (8) . The halo 
mass function is then obtained by substituting this expres- 
sion for J" into equation (9). 

2.3 Stochastic barrier 

The constant barrier Sc — 1.686 in the spherical collapse 
model is a significant oversimplification of the complex dy- 
namics leading to halo formation and growth. Such a model 
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can be improved in various ways. For instance, the excursion 
set theory results for the mass function have been shown to 
match more closely those from A^'-body simulations by con- 
sidering a moving barrier whose shape is motivated by the 
ellipsoidal collapse model (Sheth & Tormen 1999; Sheth et 
al 2001; Sheth & Tormen 2002; De Simone et al. 2010). The 
equations are summarized in Table 1. The parameters a, b, c 
are fixed by fit to A'^-body simulations, while A is fixed by 
the normalization condition on the halo mass function. As 
already remarked in Sheth et al (2001), these expressions can 
be obtained from a barrier shape that is virtually identical 
to the ellipsoidal collapse barrier, except for the factor of a, 
which is not a consequence of the ellipsoidal collapse model. 
In fact, the ellipsoidal collapse model reduces to the spheri- 
cal collapse model in the large mass limit. In this limit the 
mass function is determined by the slope of the exponential 
factor, so even in an ellipsoidal collapse model we must have 
a = 1, as in the spherical model. However, numerical simu- 
lations show that a < 1 and its precise value also depends 
on the details of the algorithm used for identifying halos in 
the simulation, e.g., the link length in a friends-of- friends 
(FOF) halo finder, or the critical overdensity in a spherical 
density (SO) finder. 

A physical understanding of the parameter a is given by 
a second independent improvement of the spherical collapse 
model, the diffusing barrier model proposed in Maggiore & 
Riotto (2010b). These authors suggested that at least some 
of the physical complications inherent to a realistic descrip- 
tion of halo formation, which involves a mixture of smooth 
accretion, violent encounters and fragmentations, can be in- 
cluded in the excursion set theory framework by assuming 
that the critical value for collapse is itself a stochastic vari- 
able, whose scatter reflects a number of complicated aspects 
of the underlying dynamics. In the simple example of a bar- 
rier performing a random walk with diffusion coefficient Db 
around the spherical collapse barrier, one finds indeed a mass 
function in which Sc is effectively replaced by a^'^Sc, with 
a = 1/(1 + Db), while at the same time k is replaced by an 
(see Table 1). 



3 HALO BIAS 

We now apply the technique in Section 2 to the computation 
of the halo bias, including the non-Markovian corrections 
with stochastic barriers. We sketch here the main steps of 
the computations, leaving the details to the Appendix. 



3.1 Conditional probability: the Markovian case 

To compute the bias, we need the probability of forming 
a halo of mass M, corresponding to a smoothing radius R, 
under the condition that the smoothed density contrast on a 
much larger scale Rm has a specified value Sm — S{Rm)- We 
use J-{Sn\Sm,Sm) to dcuote the conditional first-crossing 
rate. This is the rate at which trajectories first cross the 
barrier at (5 = (5c at time S„, under the condition that they 
passed through the point 5 = 5^ at an earlier time Sm- 
We also use the notation T{S„\0) = T{Sn\Sm = 0,Sm = 
0), so J^{Sn\0) is the first-crossing rate when the density 
approaches the cosmic mean value on very large scales. 

The halo overdensity in Lagrangian space is given by 



(Kaiser 1984; Efstathiou et al. 1988; Cole & Kaiser 1989; 
Mo & White 1996; see also Zentner 2007 for a review) 



(18) 



^ '^"'° " TiS„\0) • 

In a sufficiently large region, we have S^. <^ Sn = S and 
5m <^ 5c. Then, using the first crossing rate of excursion set 
theory and retaining only the term linear in 5m, we obtain 



l^halo 



1 



(19) 



where u — Sc/ct- After mapping to Eulerian space, one finds 
i5haio ~ 1 -|- 5^aio in the limit of small overdensity 5m — 5 
(Mo & White 1996), and 

2 _ 1 

bh{jy)^l + ^ . (20) 



3.2 Non-Markovian corrections 

We now use the path integral formalism discussed in Sec- 
tion 2.2 to compute the non-Markovian corrections to the 
halo bias. The relevant quantity for our purposes is the con- 
ditional probability 

P{5„,S„\5m,Sm)= (21) 

/_^ d5i- ■ ■ d5m ■ ■ ■ d5n,-xW {5o = 0; 5i, . . . , 5„; S'„) 

Z*;;^ d5i • • • d5m-iW (5o = 0; 5i, • • • , 5m; Sm) 

where the hat over d5m means that d5m must be omitted 
from the list of integration variables. The numerator is a 
sum over all trajectories that start from Jo = at S* = 0, 
have a given fixed value 5m at Sm , and a value 5„ at Sn , while 
all other points of the trajectory, 5i, . . . , 5m-i,5m+i, ■ ■ ■ 5„-i 
are integrated from — oo to 5c . The denominator gives the 
appropriate normalization to the conditional probability. 

Similarly to eq. (8), the conditional first-crossing rate 
J^{Sn\5m,Sm) is obtained from the conditional probability 
P{5„,S„\5m,Sm) using 

nSn\5m, Sm) = ' [^ d5^ dP{5r., S„\Sm , Sm) ^^2) 

In the Gaussian and Markovian case, the probability W^'^ 
satisfies 

W'''^{So■,5^,...,5„■,Sr^)=W'''^i5o■,5^,...,5m■,Sm) 

xW'"^{5m;5m+i,...,5„;S„-Sm), (23) 

and P{5n, Sn\5m, Sm) in eq. (21) becomes identical to the 
probability of arriving in 5n at time Sn, starting from 5m at 
time Sm, which is given by eq. (7) (with 5m identified with 
5o ai S — S„ — Sm), and we therefore recover the excursion 
set theory result. 

We have computed the non-Markovian corrections to 
this result for the case of Gaussian fiuctuations and a tophat 
filter in coordinate space, i.e. with the two-point function 
given in eqs. (15) and (16). The computation is quite in- 
volved, and we leave the details to the Appendix. Taking 
finally Sm ~ and developing to first order in 5m = 5o, 
which is the case relevant to the computation of the bias, 
for the conditional first crossing rate we find 



T{S\5o,So = 0)^ 



-Si/(2S) 



2tT 53/2 
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Table 1. Summary of Mass Function and Halo Bias Predicted by Various Analytic Models 



Model 



Mass Function f{u) 



Halo Bias bh(y) 



Parameters 



Spherical Collapse 
Ellipsoidal Collapse 

Non-Markovian 



Non-Markovian 

-I- Stochastic Barrier 



V exp I — 



'c^pi~^\ [l + (au^)'i 



i (l-«)^exp -^ +.fr 0,i^ 



K — > an, V — > \fav 



14-' 



5c = 1.686 



1 -I- -^ [v^(ai/2) -I- v^fe(azy2)i-<= A = 0.322, q = -0.3 

a = 0.707, b = 0.5, c = 0.6 



(£±11! 



"(ai/^)<: + i,(l-c)(l-c/2) 



S^ I 1-K+f e''^/2r(0,i/2/2) I 

ff [2-exp(i^)r(o,4)]} 



•{ j^^ — 1 K = for tophat-fc filter 



K = 0.35 for Gaussian 
K = 0.44 for tophat-x 



K — ^ aK, V — > \fa,v, 5c — > \/a5c a = , , „ — 

V ' V l + i-'B 

Da = diffusion coefficient 



{[l-.+ ^e^'^^muV2)) 



(24) 



5c 



K-l)+2(2-e 



X0,^V2))]} 



From this we obtain the Lagrangian halo bias b^ and the 
Eulerian halo bias bh'- 



bh{u) 



l + &h 



(25) 



1 + 



5c l-K+fe'''/2r(0,z/V2) 



{. 



1/^/2 



r(o,i/72) 



For K = we recover the usual Markovian result in equa- 
tion (20). For V ^ 1, corresponding to large masses, 
e" "r(0, !^^/2) — > 2/z^^, and the above expression simpli- 
fies to 



6h(t 



1 + 



i p 

1 — K (5c 



1-K+ 4- 
for 1/ > 1 



(26) 



3.3 Adding a stochastic barrier 

In the presence of the stochastic barrier described in Sec- 
tion 2.3, we can easily modify the halo bias in equation (25) 
using the substitution (5c — >■ a^ (5c and k — > ok, where the 
parameter a is related to the diffusion coefficient of the bar- 
rier. The Eulerian halo bias finally reads 



hh(y) 



1 + 



1 



V^(5c 1 - aK + i^e'^'''/2r(0, 01/2/2) 



(27) 



|(ai/ - 1) + Y 



=""'/^r(0, ai.V2) I 



We note that equation (27) raises the halo bias in the large 
V (i.e. large halo mass) region compared to the bias in the 
ellipsoidal collapse model (Sheth et al 2001), and in fact get 
closer to the spherical result. For i/ ^ 1, using again the 
asymptotic expression of the incomplete Gamma function 
and keeping only the leading term ~ v^ , our result reads 



hh{p) 



(1 — an) &c 



for !/ > 1 



(28) 



which differs from the asymptotic spherical collapse result 
by an overall factor of a}^^ jiX — uk). 



4 COMPARISONS 

We now compare the predictions for the Eulerian halo bias 
&h(i/) from our non-Markovian and stochastic barrier model 
with those from the standard excursion set theory as well 
as A''-body simulations. We present the results for the halo 
mass function in parallel since as we have shown in Sec- 
tions 2 and 3, an analytic theory for halo formation provides 
simultaneous predictions for the mass function and bias. 

Our model contains two parameters: (1) k, which pa- 
rameterizes the degree of non-Markovianity and its exact 
value depends on the filter function used to smooth the den- 
sity, e.g., K = 0,0.35,0.45 for tophat in momentum-space, 
Gaussian, and tophat in coordinate filter, respectively; (2) 
a, which parameterizes the stochasticity of the diffusing bar- 
rier with the diffusion coefficient Db, where a = 1/(1-1- Ds). 
There is no a priori reason to favor one filter to another, 
nor is the choice of filters limited to the three functional 
forms given above. Furthermore, we recall that in Maggiore 
& Riotto (2010a) the scaling k — >■ a k is obtained under the 
simplified assumption that the barrier makes a simple Brow- 
nian motion around the spherical collapse barrier; for more 
complicated stochastic motions of the barrier (and also for 
fiuctuations around the ellipsoidal barrier), the rescaling of k 
might be different. For these reasons, we prefer in this work 
to treat both a and k as free parameters and use simulations 
to calibrate their values. 

For the simulations, we choose to compare with the fits 
to the latest iV-body simulations (Tinker et al 2008, 2010). 
These papers provide detailed discussions about the com- 
parison within the different A^-body results and numerical 
issues such as the dependence of the results on simulation 
resolution and halo definitions and finders. 

Fig. 1 shows the results for the bias (left panels) and 
mass function (right panels) as a function of halo mass (as 
parameterized by i/ = 5c/ cr). The bottom panels show the 
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Figure 1. Comparison of tlie Eulerian lialo bias hy^iy) (left panels) and tlie halo mass function f{v) (right panels) from various analytic 
models and simulations: our non-Markovian and stochastic barrier model from eq. (27), with a = 0.818, k = 0.23 (black solid line), 
the standard Markovian spherical collapse (blue dotted) and ellipsoidal collapse (magenta dot-dashed) models, and the fits to Af-body 
simulation results (red dashed) from Tinker et al. (2008,2010). The bottom panels show the fractional difference between each of the 
analytic model prediction and the fit to Af-body result. With only two free parameters, our new model is able to match the Af-body 
results to within ~ 20% for a wide range of u. 



fractional difference between each model prediction and the 
fit to A'^-body results. For the spherical and ellipsoidal mod- 
els, we use the standard parameters listed in Table 1. For our 
model, we plot the predictions using a = 0.818 and k = 0.23, 
which provide a good match (within ~ 20%) to both the 
bias and mass function from A''-body, and in particular to 
the mass function at high mass. 

We note that since k is related to the two-point corre- 
lation function of the density field (see eq. 91 of Maggiore 
& Riotto (2010a)), it can in principle depend on the cos- 
mological model. As discussed in footnote 10 of Maggiore 
& Riotto (2010a), however, the dependence of k on the cos- 



mological parameters is extremely weak. Both k and a can 
therefore be treated as universal parameters whose values 
can be calibrated with simulations. 



5 CONCLUSIONS 

We derived an analytic expression (eq. 27 with a — 1) for 
the halo bias in the non-Markovian extension of the excur- 
sion set theory. This new model is based on a path inte- 
gral formulation introduced in Maggiore & Riotto (2010a), 
which provides an analytic framework for handling the non- 
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Markovian nature of the random walk and for calculating 
perturbatively the non-Markovian corrections to the stan- 
dard version of the excursion set theory. The degree of non- 
Markovianity in our theory is parameterized by a single vari- 
able, K, whose exact value depends on the shape of the filter 
function used to smooth the density field, e.g., k = for a 
tophat filter in momentum space, k « 0.35 for a Gaussian 
filter, and k ~ 0.44 for a tophat filter in coordinate space. 

As already discussed in Bond et al. (1991); Robertson 
et al. (2009); Maggiore & Riotto (2010a), changing the fil- 
ter function in the spherical collapse model from a tophat 
in momentum space to a tophat in coordinate space does 
not help alleviate the discrepancy in the mass function be- 
tween the Press-Schechter model and A*'-body simulations. 
In another word, had we plotted the corresponding curves 
in Fig. 1 using a — 1 (i.e. a constant barrier height as in 
the spherical collapse model) and k = 0.44 (for a tophat 
filter in coordinate space), the bias would be too high by 
up to ~ 80% at large u compared to the A^-body result, 
and the mass function would be too low by up to ~ 80% 
at large mass (see also Fig. 9 of Maggiore & Riotto 2010a). 
Using a Gaussian filter reduces k by only ~ 20% and has 
only a minor effect. Additional modifications to the theory 
beyond including non-Markovian corrections must therefore 
be introduced to match the A-body results. 

We have explored one such modification by allowing 
the barrier height itself to be a stochastic variable (Sec. 2.3 
and Maggiore & Riotto 2010b). This new ingredient intro- 
duces a second parameter a in our theory, as summarized 
in Table 1. As the solid black curves in Fig. 1 illustrate, 
an appropriate choice of these two parameters for the non- 
Markovian correction and stochastic barrier (e.g., k = 0.23 
and a = 0.818) produces a good match to A-body results 
for both the halo mass function and bias, with fractional de- 
viations being ~ 20% or less. In comparison, the ellipsoidal 
collapse model contains four fitting parameters (a, b, c and q; 
see Table 1) and does a comparable job at matching A-body 
simulations (dot-dashed magenta curves in Fig. 1). 

Further improvement to the model presented in this pa- 
per can be obtained by computing the bias through the ex- 
cursion set theory starting from the ellipsoidal model and 
including the effects of non-Markovianity. A step towards 
this computation has been taken recently in De Simone et 
al. (2010), where the first crossing rate has been computed 
using the path integral method for a generic barrier. As the 
next step, one could envisage to combine the diffusing bar- 
rier model with the ellipsoidal model, i.e. consider a barrier 
that fiuctuates around an average value given by the ellip- 
soidal collapse model, rather than around the constant value 
provided by the spherical collapse model as done in this pa- 
per. We expect this combined model to be able to provide 
an even closer match to A-body results than the ~ 20% 
accuracy achieved by either model alone. 
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APPENDIX A: DETAILS OF THE 
COMPUTATION 

To perform the computation we use the technique discussed 
in detail in Maggiore & Riotto 2010a (MRl). We consider 
first the numerator in eq. (21). We start from eq. (12), with 
the two-point function {5i5j)c given in eqs. (15) and (16), 
and we expand to first order in n (recall that Aij is propor- 
tional to k). This gives W in terms of W^'^, 

W(5o;...,5n;S„) = 



~W^'"{So;...,S„;S„) 

n 

+ i Y^ A,,9,9,iy«^-"(5o; • • • , 5n; S^) , 



(Al) 



j,i=i 



where A; 
identity 

AfeC 



A{Si,Sj), di = d/dSi, and we have used the 



EU^^'^ 



= -idkC 



EU'^' 



(A2) 



to transform the factor —AijXiXj coming from the expansion 
of the exponential into Aijdidj. It is convenient to split the 
sum into various pieces 



n m— i 771—1 

- ^ A^jdidj = - ^ Aijdidj + ^ A 



did„ 



im'^i^m 



»J = 1 



*.i=i 
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ri — ± ri. — 1 

i,j— rri + 1 i— rri + 1 



m—1 n—1 



+ 2J 2_/ ^^J'^i'^3 + 2_/ ^3"zdjd„ 



i — 1 J — 771 + 1 



J — m. + 1 



(A3) 



Consider first the contribution from the first line of this 
expression. Using the factorization property (23) of 14^^™, 
its contribution to the numerator in eq. (21) can be written 
as 

fSa 

dSi ■ ■ ■ dSm-idSm+i ■ ■ ■ dS„-i 






i,j^l i—l 



xW'^iSo; ..., Sm;Sm)W<''^(5m; . . . , cS„; S„ - Sm) 



m—1 m— L 



m—1 777 — 1 

= / dSi- ■ ■ dSm-1 

X / d5m + l ■ ■ ■ dSr.-lW'^iSn,; . . . , <5„ ; Sn - S™) 
•^ — oc 

dSl--- dSm-l Y^ /:^^ra^^W^'^{So■, . . . , Srr.] Srr.) 

/.5c 
dSm + 1 ■ ■ ■ d5n-ldmW^'^ {5m; . . . , S„; S„ - Sm) 
-oc 



(A4) 



The first term is easily dealt by observing that 

dSm+i ■ ■ ■ dS„-iW^"'{Sm; . . . , (5n; 5„ — Sm) 



(A5) 



Combining this with the contribution coming from the zero- 
th order term W^'^{5o; ■ ■ ■ ,Sn; Sn) in eq. (Al) and using 
again the factorization property (23) of W^™, we therefore 
get 

Il^'^{5m;5n;Sn — Sm) / dSl ■ ■ ■ dSm-l 

J —oc 

m — 1 



1 + - 2^ Aijdidj + 2^ Aimdidn 



».i=i 



xW'''^i5o;...,5m;Sm) 



(A6) 



+ / d5i---d5m-i^A,mftW'^'"(5o;...,c5m;5'm) 
X / dSm+i---dS„^idmW^"'{Sm;- ■ ■ ,Sn;Sn- Sm) ■ 



We now observe that the term in brackets give just the ex- 
pansion to 0{k) of the denominator in eq. (21). Therefore, 



to 0{k), we can write 

^\0n^ On\0m, Om) -— 11 yOm] On] On Om) 



+ P 



non — mark 



{Sn,Sn\Sm,Sm) , (A7) 



where 

P {Sn,S„\5m,Sm)= TT^m ^ f : ET^ . (^8) 



ns'^{5o;Sm;Sm) 



and Na , ■ ■ ■ , Nd are defined by 

fSa m-l 



Na = 



d<5l ■ ■ ■ dSm-l J2 A^m^,W<''^{So■, ...,Sm;S„ 

'OO ■ 1 



X / dSm + l-'-dSr^-ldmlV'^^iS. 



^n 1 t-'77, ^rr 



TTi ) • • • 1 '-'77 1 ^-^n ^m) 5 



N. = 



,5c 



dSi ■ ■ ■ dSm-idSm+i ■ ■ ■ dSn-1 



(A9) 
(AlO) 



- 2J Aijdidj + 22 ^^ridtdn 

7,J— TTl+l i — TTl+l 

X VF'^'"(r5o; . . . , 5m; Sm)W<''^{5m; . . . , r5„; S„ - Sm) , 



Nc = j dSl- ■ ■ dSm-ldSm+l ■ ■ ■ dSn-l (All) 

— 00 
m—1 

xW'^iSo; ...,Sm; Sm)W'"^{5m; . . . ,5„;S„ ~ Sm) , 



Nd = 



d^i • • • dSm-idSm+i ■ ■ ■ dS„-i 



(A12) 



777—1 n — 1 



22 z2 ^»J^»^J + 2_/ ^]rndjd„ 

7^1 J— TTl+l J— TTl+l 

xW'^iSo; ■■■,Sm; Sm)W^'^{5m; ...,5n;S„- Sm) • 

The contribution Na comes from eq. (A4), while Nb, Nc and 
Nd come from the second, third and fourth line in eq. (A3), 
respectively. Observe that in the denominator in eq. (A8) we 
could replace n{5o;Sm; Sm) by n^™(5o; (5™; 5*™), since the 
numerator is proportional to Aij and therefore is already 
0{k). 

The contributions Na , ■ ■ ■ , Nd can be computed using 
the techniques developed in MRl. The term Na is imme- 
diately obtained using eqs. (105) and (110) of MRl, and is 
given by 

^ J.i5c-Sm) ^^^j25c-5„ 



•Jm \ sj ZOrr 

X Omii yyra '•, On \ On *-'777 J • 



(A13) 



where Erfc is the complementary error function. The term 
Nb is given by 

= W-^{kv,&m;Sm) (A14) 

X[n (<5m, S'm;<5n,5'n) + n {&m, Sm\&n, Sn)\ , 



N, 



where 

11 [Om^ Om.; On, On) 



dSm + l ■ ■ ■ dSn-l (A15) 

CO 
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n-l 

Y^ A,„d,d„W^'^(Sm;...,S„;S„-Sm), 



i — TTl + l 



and 



(A16) 



n {Sm,Sm.;S„,S„) = / dSm+1 ■ ■ ■ dSn-1 
n-l 

The computation of n''^ and H*^ is quite similar to the com- 



putation of the terms called n™"™ and W 
and in the continuum limit e — >■ we get 



mem — mem 



n {Sm,Sm;Sn,S„) ^ d„ lira - 



dS, 



X A(5'„ S„)nnSm; 5c; S, - 5„.)nr(5c; S„; S„ - S^) 

fSr, 
— {5c — 5m,)dn,l {5c — 5n) / dSi 



{5c ~ Sm) {Sc — Sn) 



X exp 



^\Oi Om) '^\^n ^i) 



(A17) 



and 



dSr / dS, 



n {5m.,Sm;5n,Sn) ^VaXl — 



xA{S„Sj)Iir{5m\5c;S, - S,^) 
xnr{5c; 5c; S, - S^nr (5c; c5„; S„ - S,) 



K 



: {5c — 5m){5c — 5n) 






s, 



-(ic-<S™)V[2(S,-S„)] 



dS, 



5,(S, -5Oi/2(5„-S,03/2 



(A18) 



This can be rewritten as a total derivative with respect to 
5„, as 






(5c — 5m)dn 



"" 'L '^^^ '(S^ - SmY''' 



7r\/27r 



dS,- 



g-(<S<:-5„) = /[2{S„-S,)] 



(A19) 



The fact that both II''^ and II'''^ can be written as a deriva- 
tive with respect to 5„ simplifies considerably the compu- 
tation of the fiux J-{S), since we can integrate dn = d/d5n 
by parts, and then we only need to evaluate the integrals in 
eqs. (A17) and (A19) in 5„ — 5c, which can be done analyt- 
ically, as discussed in MRl. 

The term A^c is a total derivative with respect to d„ of 
a quantity that vanishes in 5n ~ 5c so, when inserted into 
eq. (22), it gives a vanishing contribution to the first crossing 
rate. The most complicated term is Nd- Using the techniques 



developed in MRl, a rather long computation gives 

25c — Sn 



Nd = -dJ5c{5c - 5^)Eric ( ^^^^-^ ) 



dm I 



(A20) 
(A21) 



in MRl, J(5„,5„) 



V2S„ 
+ n<"^{5iy,Sm;Sm)Ndy 
where 

Ad = —5m{5c — 5m)I{5m,5n) + Sm{5c — 5m)dml{5m,5„) 
~SmI{5m,5n), (A22) 

and 



dS". 



X exp 



's,{s,^s^y/^Sr^~s,Y/^ 

{5c - 5mf {5c - 5r.f 



(A23) 



^(Oj Orn) '^y^n ^j) 

Using eqs. (A8) and (A14), cq. (A7) can be rewritten as 

P(5„, s„\5m, Sm) = n^-+n"+n''^+ J^:+^; + f .(A24) 

lls™(do;dm;6m) 

We can now compute the contribution to the flux from the 
various terms. The term H^™ gives the zero-th order term, 

Obn ' 



Sc — Sn 



' — oc 
,-(6a-Sm.f/l2(Sr,-Sm.)] 



The contribution of n''^ to the fiux is zero since it is the 
derivative with respect to 9„ of a quantity that vanishes in 



•^ yJnl'^m, *Jm.} 

d5„Il {5m,Sm;5n;Sn) 



(A26) 



d 
'dS„ 



d 



-K^/l-K dSn 



X 



dS, 



dS, 



S", 



{S^ - S„)3/2 



-(.5,-«„)2/[2(S,-S„)] 



1 



I Si Sj{bj bi) / {bn bj) 

The inner integral is elementary, 

r-S„ 1 

dS, 



g^ 5,(5, -S',)i/2(S„-S,)i/2 (5.5„)i/2 

and we end up with 



J~ ybnl^m.'i bm) 



d 

'dS„ 



X / dS^ 



5 



1/2 



ybi bm) ' 



' k{5c — 5m) 
- \/27r 5n 



(A27) 



(A28) 



(A29) 



which generalized eq. (118) of MRl to 5m 7^ and Sm 7^ 0. 
For Sm generic the integral cannot be performed analyti- 
cally. However, for computing the bias we are actually in- 
terested in the limit Sm — >■ with 5m generic, and we see 
that in this case this contribution reduces to that computed 
in MRl, with the replacement 5c ^ 5c — 5m- 
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The remaining contributions can be computed similarly. 
For the term Nd, again, rather than computing explicitly the 
derivative dn = d/dSn in eq. (A20) , it is convenient to insert 
this expression directly into the first-crossing rate (22), and 
use the fact that it is a total derivative with respect to dn to 
perform the integral over dS„. So, in the end, we only need 



I{Sm, Sn = Sc) = 



„+('5c-«™)V(2S„) 



yJm^n) ' 



(A30) 



xErfc 



(<5c — <5m) 



On 



^Oni\On ^vfi) 



We can now put together all the terms and take the limit 
Sm — > (with 5c — <5m > 0). In this limit the Erfc function 
in eq. (A30) reduces to an exponential, so 



/(5m,l5„ = (5c) 



\j2-K(Sn — Sm) 



-(.5,-«„)2/[2(S„-S„)] 



(Sc — 5n)Sn 

Denoting 5m = (5o in this limit, we finally get 
1 — K (5, 



(A31) 



T{S\5o,Sm = 0) = 



+ - 



2^ 53/2 



^g-(^c-.So)"/(2S) 



5o 



"^"r fo ('^'=^'^°)^ 



2^53/2 



2S 

(<5c-<5o)' 



(A32) 



S 



„-(ic-<So)^/(2S) 



Expanding this result to first order in 5o we obtain eq. (24). 
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